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ABSTRACT 

We present an axisymmetric, equilibrium model for late-type galaxies which 
consists of an exponential disk, a Sersic bulge, and a cuspy dark halo. The model 
is specified by a phase space distribution function which, in turn, depends on 
the integrals of motion. Bayesian statistics and the Markov Chain Monte Carlo 
method are used to tailor the model to satisfy observational data and theoret- 
ical constraints. By way of example, we construct a chain of 10 5 models for 
the Milky Way designed to fit a wide range of photometric and kinematic ob- 
servations. From this chain, we calculate the probability distribution function 
of important Galactic parameters such as the Sersic index of the bulge, the disk 
scale length, and the disk, bulge, and halo masses. We also calculate the probabil- 
ity distribution function of the local dark matter velocity dispersion and density, 
two quantities of paramount significance for terrestrial dark matter detection 
experiments. 
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Though the Milky Way models in our chain all satisfy the prescribed ob- 
servational constraints, they vary considerably in key structural parameters and 
therefore respond differently to non-axisymmetric perturbations. We simulate the 
evolution of twenty-five models which have different Toomre Q and Goldreich- 
Tremaine X parameters. Virtually all of these models form a bar, though some, 
more quickly than others. The bar pattern speeds are ~ 40 — 50kms _1 kpc~ x at 
the time when they form and then decrease, presumably due to coupling of the 
bar with the halo. Since the Galactic bar has a pattern speed ~ 50kms~ 1 kpc _1 
we conclude that it must have formed recently. 

Subject headings: Galaxy: kinematics and dynamics — methods: statistical - 
methods: N-body simulations — cosmology: dark matter 



1. INTRODUCTION 



Dynamical galactic models serve a variety of purposes. They may be used to interpret 
the structural and kinematical observables of galaxies - surface brightness profiles, rotation 
curves and velocity dispersion profiles - in terms of intrinsic three-dimensional density and 
velocity distributions. Dynamical models also provide a starting point for controlled simula- 
tions of complicated processes such as the formation of bars and spiral structure. In short, 
galactic modeling provides the essential interface between observations and detailed theories 
of galaxy formation. 

In this paper, we introduce a new dynamical model for late-type galaxies which com- 
prises a disk, bulge, and dark halo. The model is derived from equilibrium solutions to the 
collisionless Boltzmann and Poisson equations. It extremely flexible and may be tailored to 
satisfy observational data and theoretical constraints. We use Bayesian statistics and the 
Markov Chain Monte Carlo (MCMC) method to implement these constraints and to deter- 
mine the probability distribution function (PDF) of the model in the full multi-dimensional 
parameter space. 



Our model builds upon earlier work by lKuijken fc Dubinskil (119951 ) and lWidrow Dubinski 



( 120051 ). The original Kuijken & Dubinski model consists of an exponential disk, a King-model 
bulge, and a lowered Evans-model halo and has the attractive feature that the phase space 
distribution function (DF) is built from analytic functions of the integrals of motion. No addi- 
tional assumptions about th e velocity-space d istribution are made. By con t rast, t he widely- 
used approach described in iHernquistl (119931 1 (see also ISpringel fc White! (119991 )) assumes 
that the local velocity distributions of the halo and bulge particles are Gaussian with disper- 
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sions estimated from the Jeans equations. This approach leads to models which are slighly 
out of equilibrium. When used as initial conditions in N- body experiments, they readjust to 
a diffe rent state from the one proposed (see, for example, iKazantzidis. Magorrian. &: Moore 
( b004h ). 



There are two main disadvantages of the Kuijken and Dubinski models. First, the 
bulge and halo have constant density (or weakly cuspy) centers whereas actual bulges and 
dark halos may have central density cusps. Second, the structure of the bulge and halo 
are determined implicitly by the model parameters. (By construction, the disk's structural 
parameters, namely its radial and vertical scale lengths, its mass, and its truncatio n radius, 
are determined explicitly by the model parameters.) IWidrow &: Dubinski! (120051 ) built a 



galactic model with r •'■-density cusps for both the bulge and halo but again, with a DF that 
determines the structure of the bulge and halo implicitly. 

For our new model, the closed-form DF is abandoned in favor of a numerical DF which 
is designed to yield, to a very good approximation, user-specified density profiles for the 
bulge and halo. That is, the density profiles of the bulge and halo are now explicit functions 
of the model parameters. The present version of the model allows for a Sersic bulge and a 
halo with p oc r -7 as r — > where 7 is between and 2. 

Our model is specified in terms of fifteen or so parameters. How are these parameters 
selected? One approach is to choose models at random and identify the ones that satisfy 
certain general constraints (e.g., the Tully-Fisher and size-luminosity relations). The result 
would be a catalog of disk-bulge-halo systems which could be used to study kinematical and 
dynamical trends such as the circular speed-central velocity dispersion (V c — <tq) relation 
(see, for example, ICouteau et al.l (120071 )). A catalog of this type could also be used for N- 
body studies of mergers or interactions between galaxies in a cosmological environment. A 
second approach, and the o ne pu rsued here, is to build mode ls for specific galaxies. Following 
Kuijken &: Dubinski! (119951 ) and IWidrow fc Dubinskil (120051 ) , we use the Milky Way as our 
illustrative example. 



M odeling the Milky Way is a time-honoured endeavor; notable exampl e s include llnnanen 



( 1973 ) 
(1981) 



Clutton-Brock. Innanen. fc Pappl (1977 ). Bahcall fc Soneira ( 198of). Caldwell fc Ostriker 



Kuiiken fc Gilmore! Jl99lh iRohl fe fc KreitschmanrJ dl988h.lMaJhotral Jl995h:lKochanekl 

Jl996h . lEvans fe Wilkinson] Jl999h . and iKlvpin. Zhao, fc Somerville! (booi l. Our construc- 



tion o f dynamical Milky Way models is in the spirit of the mass model survey by lDehnen fc Binney 
(119981 ) . Their models comprise a multi-component disk, a bulge, and a halo and are charac- 
terized by ten parameters. Twenty-two examples are presented, each the result of a maximum 
likelihood analysis in which some parameters are held fixed while others are allowed to vary. 
The Dehnen & Binney likelihood function is constructed by comparing model predictions 
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with six sets of observational data: the inner and outer Galactic rotation curves, the Oort 
constants, the mass at large radii, the local vertical force, and the line-of-sight velocity dis- 
persion in Baade's window. The advantage of our models is that they not only describe the 
potential-density pair for the Galaxy but also, the underlying DF. We therefore have the 
ability to examine the stabili ty of our Galact ic m odels using N-body simulations, an issue 
that is often ignored (but see ISellwoodl (119851 ) and iFuxi (119971 )). 



For the most part, we adopt Dehnen & Binney's choice of observational data though we 
include more complete observations of the line-of-sight dispersion in the bulge region as well 
as photometric data from the COBE satellite. We also present what we believe to be a more 
balanced treatment of the likelihood function. Most significantly, we bring to the problem 
the powerful tools of Bayesian statistics and MCMC. These tools allow us to map out PDFs 
of both input parameters and derived quantities. 

Though our model represents an axisymmetric, equilibrium system, it is susceptible to 
non-axisymmetric instabilities and therefore provides a natural starting point for numerical 
studies of galactic dynamics. An N-body realization of the model can be easily generated 
from the DF and then used as the initial conditions for a numerical simulation. 

The Milky Way models in our MCMC series all satisfy the observational constraints 
but vary considerably in their structural propert ies. We simul ate a selection of twenty-five 
models which span a wide range i n Toomre Q ( IToomreill964l ) and Goldreich-Tremaine X 
( IGoldreich &: Tremaindll978l . Il979l ) parameters and find that a bar develops in virtually all 
of the cases. The onset of the bar instability can occur immediately or after several Gyr, 
depending on the model. 

We present the model in Section 2, review the observational constraints in Section 3, and 
provide a summary of the essentials of Bayesian statistics and the MCMC method in Section 
4. We discuss some preliminaries including our choice of prior probabilities, in Section 5. We 
present the results of our MCMC analysis in Section 6 and the results of our bar formation 
simulations in Section 7. In Section 8 we summarize our main conclusions and speculate on 
how we might improve upon and extend the models and MCMC analysis. 



2. GALACTIC MODELS 



We consider axisymmetric, collisionless systems whose DF is of the form 

f(£,L z ,E z ) = /disk (£, L z , E x ) + /bulge (£) + /halo (£) 



(1) 
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where £ = —E is the relative energy, L z is the angular momentum abou t the symmetry axis, 
and E z is the energy associated w ith vertical motions of stars in the disk (IKuijken Dubinski 
19951 ; IWidrow fc Dubinskill2005l ). For time- independent, axisymmetric systems £ and L z are 
integrals of motion while E z is an approximate integral of motion for disk stars on nearly 
circular orbits. Jeans theorem implies that a system generated by equation [1] will be in 
approximate equilibrium. 

Integrating equation^] over all velocities yields the density in terms of the gravitational 
potential, $, and the cylindrical coordinates R and z: 



P {R-i z i ^) — Pdisk (R, z i ^) + Pbulge (^) + Phalo 



(2) 



where $ = — $ is the relative potential. (Note that implicit in equation [T] is the assumption 
that the bulge and halo velocity dispersions are isotropic.) Self-consistency requires that p 
and \l/ satisfy Poisson's equation 



V 2 * 



A-np (R, z, ^) 



(3) 



which is accomplished, in practice, through an iterative scheme. (Note that here and through- 
out, we set Newton's constant (7=1.) 



Kuijken &: Dubinski ( 1995 ) chos e fhnige to be the King model DF ( Kingj[l966 ) and /halo to 



be the lowered Evans model DF from lKuijken &: Dubinski! (I1994I ). (The latter depends on L z 
as well as £ thereby allowing for flattened halos.) Their models have two main shortcomings. 
First, bulges and halos may have central density cusps whereas the King and lowered- Evans 
DFs yield density profiles with constant-density cores. Second, the relationship between 
the model parameters and the density profiles of the bulge and halo is implicit rather than 
explicit and not particularly intuitive. 



Widrow fc Dubinski! (120051 ) built g alactic models wi th cuspy (p oc r 1 as r — > 0) bulges 



and halos. Sp ecifically, they chose the iHernquistl (119901 ) DF for the bulge and a DF from 
Widrowl (120001 ) for the halo. The latter was constructed to yield the so-called NFW profile 



PNFw( r ) = — 2 (4) 

(r/a h ) (1 + r/a h ) 

rtNavarro. Frenk fc Whit3l996h . I Widrow fc Dubinskil (120051 ) altered these DFs in an attempt 
to account for the gravitational potential of the other components though the modifications 
were somewhat ad hoc. Consequently, the density profile of the halo differed from the NFW 
profile and likewise for the bulge. 



A further drawback of the lWidrow fc Dubinskil (120051 ) model is that the Hernquist bulge 
and NFW halo are arguably too restrictive. The Hernquist DF yields a system whose surface 
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density profile is well approximated by the i? 1 / 4 -law ( Ide Vaucouleurall948l ). The bulges of 
late-type galaxies are found to have surface brightness profiles which follow the more general 
Sersic law, 

S(r) = S oe - fc(W)1/ " , (5) 



with Sersic index n between 0.6 and 2 (lAndredakic. Peletier. fc Balcellslll995l ; ICourteau. de Jong, fc Broeils 
19961 ). Likewise, dark matter halos may have density p rofiles more general than the NFW 
form. Since the work oflNavarro. Frenk fc White (1996). there has been considerable debate 



Moore et ali (19991) find e viden ce in their sim- 



over the actual form of halo density profiles, 
ulations for steeper cusps (p oc r~ L5 ). More recently, iNavarro et al.l (J2004]) conclude that 
the logarithmic slope of the halo density profiles decreases steadily with radius though their 
results are still consistent with equation HI On the observational front, the rotation curves 
of dar k matter-dominated low surface brightness galaxies appear to favour constant density 



cores (lMoorelll994l ; iFlores fc Primacklll994l ; iMcGaugh fc de Bloklll998l ; Ivan den Bosch et al. 



20001 ) though this interpretation of the data is being challenged on a number of fronts. 



For our new models, we begin by choosing target density profiles, pbuige and phaio, for 
the bulge and halo. Assume, for the moment, that the system is spherically symmetric. 
Through the Abel integral transform, 



l 



8tt 2 



d 2 p 



total 



d*tatal V£ ~ * total 



bulge or halo 



(6) 



(IBinney fc Tremaind 119871 ). we can construct bulge and halo DFs which yield the target 
density profiles. In the case of an isolated halo or bulge, ^totai is the potential derived from 
Phaio or pbuige and equation[6] reduces to the usual expression for the DF of a spherically 
symmetric system with isotropic veloci ti es. T he DF for a system following the Sersi c law 
was found with this m etho d by ICiottil (19911 ). DF s for NFW halos were found by IZhao 
rtl997h . lWidrow1 (boooh . and lLokas fc Mamonl (120001 ) . For a composite system or one with 
an externa l potential, one s i mply replaces the $ derived from pi with the total gravitational 
potential. iTremaine et al.l (120021 ) used this method to derive DFs for bulges with central 
black holes by setting * total = ^ bulge + GM b i ackho i e /r. 

Equation[6] is only valid for spherically symmetric systems, a condition violated once a 
disk is included. Our approach is to use a spherical approximation (essentially, the monopole 
term of a spherical harmonics expansion) for the disk potential in evaluating \1/ total . We stress 
that equationElis used to construct /haio(£) and /bui ge (£), not to solve for ^(R, z) and p(R, z). 
We can use /haio(£) and /bulge (£) in equation [1] even though the composite system is not 
spherically symmetric; A DF of the form / = /(£) yields an equilibrium system in any time- 
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independent potential regardless of the spatial symmetries of the potentialjj The bulge and 
halo of the final model are axisymmetric (but not spherically symmetric) since isodensity 
surfaces follow isopotential surfaces (equation |2j), the latter being flattened by the disk. As 
we demonstrate below, the spherically-averaged density profiles of the bulge and halo are 
very close to the target profiles. 



2.1. Target Density Profiles 
We choose the target density profile for the bulge to be 



-p 



/Wr) = *(JJ e -*/*-> l/ " (7) 

which yields the Sersic law (equation [5J for the projected surfac e density profile provide d 
one sets p = 1 - 0.6097/n + 0.05563/n 2 JPrugniel fc Simienl \l99^ iTerzic fc Graham! l2005h . 
Note that here and in equation [5], Eq, R e and n are free parameters while the constant b is 
adjusted so that R e encloses half the total projected light or mass. In our models, we use 

a b = (A-Knb n ^T (n (2 - p)) R 2 e p b ) 1/2 (8) 

rather than p& to parametrize the overall density scale of the bulge models. With this 
definition, <j\ corresponds to the depth of the gravitational potential associated with the 
bulge. 

We choose the target density profile of the halo to be 

2 2 ~ 7 ct 2 1 

Phalo = , 2 7 3-7 C ( r ? T ^ Sr h) (9) 

4vra h (r/a h ) (1 + r/a h ) 

where C is a truncation function that smoothly goes from unity to zero at r = r% over a 
width 5rh- We use the function 

c < ri ^H erfc te) ■ (10) 

The halo profile is therefore characterized by five parameters: r^, 5rn, the halo scale length, 
ah, the velocity scale, cr^, and the central "cusp strength", 7. For r < r h , the mass interior 
to radius r is given by 



1 A self- consistent syste m with a DF that depends only on the energy must be spherically symmetric 
( Binnev fc Tremainel 1987 ). The statement in the text applies to systems in which there is an external 
potential that does not necessarily respect spherical symmetry. Here the potential due to the disk plays the 
role of an external potential to the halo and bulge. 
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M(r) 



I 'a h a h 



l + r/a h 



+ log(l + r/a h ) 



Following iKuij ken &: Dubinskil (119951 ). we adjust the disk's DF so that its space density 
falls off approximately exponentially in R and as sech 2 in z with radial and vertical scale 
lengths Rd and Zd respectively. The disk is truncated at a radius R ont with a truncation 
sharpness parameter SRd- In addition, we choose a DF where the radial dispersion profile is 
approximately exponential: 



o 2 R {R)=o* a exp(-R/R ir ) 



(12) 



For simplicity, we set R a = Rd in accord with observations by iBottemal (119931 ) 



OBSERVATIONAL CONSTRAINTS 



We use nine sets of observational data as constraints on our Milky Way models. Five 
of these data sets - the inner and outer rotation curves, the Oort constants, the vertical 
force in the solar n eighb orhood, and the total mass at large radii - are largely taken from 
Dehnen fc Binneyl (119981 ) and references therein. We inco rporate measurements of the line- 
of-sight bulge dispersion from the compilatio n of data by iTremaine et al.l (120021 ) as well as 



estimates of the local velocity ellipsoid from lBinney fc Merrifie 



d (1998 ). We also use dust- 



corre cted near-infrared DIRBE data from the COBE satellite (IBinney. Gerhard, fc Spergel 
mi- 



Inner rotation curve 

Inside the solar circle, the Galactic rotation curve is usually presented in terms of the 
"terminal velocity" , the peak velocity along a given line-of-sight at Galactic coordinates 
6 = and |/| < ir/2. Assuming that the Galaxy is axisymmetric, v term is given by 



Vte 



v c (R) - v c (Rq) sin/ 



(13) 



is the distance from the Sun to the Galactic center an d v r is the circula r 



where R 

speed (see, for example, IBinney fc Merrifieldl ( 1998)). Following iDehnen Sz Binney 



( 119981 ) we use data from HI emission observations by 



Malhotral ( 119951 ) restricted to the 



range sin / > 0.3 so as to avoid distortions from the bar. 



• Outer rotation curve 
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The radial velocity of an object relative to the local standard of rest, v\ sv , is related to 
the circular rotation curve through the equation 



Vis 



~R, 



v c (R) — v c (Rq) ) cos b sin I 



(14) 



where R = (d 2 cos 2 b + R 2 , — 2R d cos b sin Z) 1 ^ 2 , (/, b) are the Galactic coordinates, and 
d is the distance to the object. Measurements of tv,; and di are compared to model 
estimates for W(R) and d(R) where W(R) = (Rq/ R) v c (R) — v c (Ro) = vi sr / cos b sin/. 
R is regarded as a free parameter which is adjusted to minimize 



f W(R)-Wj 



f g(g) 



— <i, 



(15) 



where Wi = v\ SVti / cos b sin Z. In what follows we use data from iBrand fc Blita (119931 ) 
with the same restrictions as in iDehnen &: Binneyl (119981 ) (i.e., I < 155° or / > 205°, 
d > 1 kpc, and < 0) so as to avoid contamination from non-circular motions. 

Local circular speed 

A further constraint from the rotation curve of the Galaxy comes from esti mates of the 
circula r speed at the solar radius, v c (R ). Here, we adopt the estimate of iReid et al. 



( 119991 ) who carried out VLB A observations of Sgr A*. Under the assumption that Sgr 
A* is at the center of the Galaxy, they found 



v c (Rq) = (219 ±20 km s" 1 ) 



Rq 
8 kpc 



(16) 



which is consistent with other other estimates (jSackettlll997l ). 



Vertical force above the disk 



Kuijken fc Gilmord (119911 ) use K dwarf stars as tracers of the gravitational potential 
above the Galactic plane thereby placing a constraint on the total surface density in 
the solar neighborhood. They find 



1^(1.1 kpc) 
2txG 



71 ±6 M pc 



-2 



(17) 



independent of the relative contributions of the disk and halo. By including additional 
constraints on the local circular speed, Galactocentric distance of the Sun, and Oort 
constants, one can ferret out t he separate contributions of the disk and halo to the 
local surface density. Doing so, iKuijken fc Gilmord (I199ll ) found 
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- J disk 



48 ± 9 M pc" 



in excellent agreement with estimates of known matter in the solar neighborhood. We 
adopt equation ITTl as the constraint on the vertical force at (R, z) = (Rq, 1.1 kpc) and 
equation [TS] the constraint on the surface density of the disk at R = Rq. 



Oort constants 
The Oort constants, 



2\R 



dvc 
OR 



and 



B = 



r=r 



R OR 



(19) 



R=R 



measure, respectively, the local shear and vorticity in the Galactic disk. Here we adopt 
the constraints 



A = 14.8 ±0.8 km s" 1 kpc 



-i 



B = 12.4 ± 0.6 km s" 1 kpc 



-i 



(20) 



from the iFeast &: Whitelockl (119971 ) analysis of Cepheid proper motion measurements 
by the Hipparcos satellite. 

Local velocity ellipsoid 

The kinematics of stars in the solar neighborhood provide important constraints on 
the structure of the Milky Way. The observation that v % ^ vl alread y tells us that the 
disk DF necessarily involves a third integral of motion ( Binnev 1987h . Our cons t raints 



for the local velocity ellipsoid ar e taken from Table 10.4 o f 
which in turn were derived from lEdvardsson et al.l (119931 ) . 
give separate values for the thin and thick dis ks. Since our models assum e a single disk 



Binnev k. Merrifieldj (Il99fih 



Binnev fc Merrifieldl (ll998f ) 



component we use a mass- weighted average (IWidrow fc Dubinskil 120051 ). 



Bulge dispersion 

Observations of the line-of-sight velocity dispersion in the direction of the bulge provide 
important constraints on the bulge parameter and, to a lesser extent, the parameters 
of the other components. We use measure ments of the line-of-sigh t dispersion between 
4pc and 1300 pc from the compilation by iTremaine et al.l (120021 ) . Since the bulge is 
triaxial, the measured line-of-sight dispersion depends on the observer's orientation to 
its principal axes. Our line-of-sight to the Galactic center is approximately 20° from the 
long axis of the bulge (IBinney. Gerhard, fc Spergellll997l ) and therefore the measured 
line-of-sight dispersion will be system atically higher than the values one would obtain 
assuming a spherical bulge. Following ITremaine et al.l (120021 ). we adjust the measured 
dispersions downward by a factor 1.07 to account for this effect. 
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Mass at large radii 

The observed velocity distribution of the Milky Way satellite system and the dynamics 
of the Magellanic Stream, together with measurements of the high- velocity tail of the lo- 
cal stellar velocity distribution, provide constraints o n the large-scale mass distribution 
of the Galactic halo. Fol l owing IDehnen Sz Binnev w ho ba se their arguments 

on analyses by iKochanek! (119961 ) and iLin. Jones, fc Klemolal (119951 ) , we adopt 



M (r < 100 kpc) = (7 ± 2.5) x 10 n M 
as a constraint on the mass distribution at large radii. 



(21) 



Surface photometry 

The distribution of stars in the Milky Way is most easily determined from observations 
in the near infrared where stellar emission dominates over dust emission. Though dust 
is more transparent at these wavelengths than in the optical, extinction due to dust is 
still significant toward the Galactic center. ISpergel. Malhotra. &: Blitzl (119961 ) produced 
extinction-corrected maps of the inn er Galaxy based on t he DIRBE data set and a 
three- dimensional dust model (see also lFreudenreichl (119981 )). iBinney. Gerhard, fc Spergel 
( 119971 ) constructed three-dimensional models for the light distribution of the disk and 
bulge based on these maps with the aim of constraining the structural parameters of 
the Galactic bar. 

Our initial goal is to construct axisymmetric models for the Galaxy. Toward this 
end, we use the surface brightness as a fun c tion o f I at mid Galactic-latitudes (3° < 
\b\ < 4.5°) from lBinney. Gerhard, fc Spergell (119971 ) (their Figure 2, lower panel) where 
the effects of the bar are not so pronounced (i.e., where their axisymmetric model 
adequately reproduces the observed surface brightness profile). 



The mass model survey of IDehnen fc Binneyl (119981 ) employs a maximum likelihood 
analysis where the likelihood function is exp (— Xdb) with 



Xdb 



w, 



N. 



out 2 
Xout 



N, 



other 2 

Xother " 



(22) 



other 



The subscripts "in" , "out" , and "other" refer to the inner and outer rotation curve constraints 
and the other constraints (e.g., Oort constants, vertical force) respective ly. The Nj are the 
numb ers of data points actually used while Wi are weights introduced by IDehnen &: Binney 
( 119981 ) to account for "the number of really independent constraints". That is, the Wi 
are meant to compensate for the fact that a quantity such as the Oort constant A has been 
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obtained from a large number of data points. iDehnen fc Binneyl (119981 ) choose W m = W ou t 
W ot her = 6 though admit that the choice of WiS are "subject to ones prejudices". 



In our view, the likelihood function should be exp (— x^ ot /2) where 



Xtot 



2,2 i_ 2 
Xin ~r Xout ' Xother 



(23) 



The fact that the Oort constant constraints are obtained from a large number of (raw) data 
points is already accounted for by the small quoted errors. Dividing xf n by N in unfairly 
short-changes the rotation curve data in favour of the Oort constant constraints, and so 
forth. 



To survey the model parameter space, IDehnen fc Binneyl (119981 ) adopt the following 
approach; fix certain parameters and maximize the likelihood function by allowing the 
remaining parameters to vary. The procedure is then repeated with the fixed parame- 
ters set to different values, or different subsets of parameters held fixed. The result is a 
rather uneven survey of t he ful l parameter space. A similar exercise was carried out by 
Widrow. Perrett. fc Suyul (120031 ) for M31 using the original Kuijken and Dubinski models 
together with rotation curve, velocity dispersion, and surface brightness data. T his proce- 
dure was also used for both M31 and the Milky Way in I Widrow fc Dubinski! (120051 ). Bayesian 
statistics and MCMC provide a more complete picture of the model parameter space as we 
now demonstrate. 



4. BAYESIAN ANALYSIS AND MCMC 

Our aim is to calculate the posterior probability density function, p(M\D, J), of a Galac- 
tic model M given data D and prior information I. M is specified by the fifteen model 
parameters as well as additional astronomical parameters - here i?o and the mass-to-light 
ratios of the disk and bulge, (M/L) d and (M/L) b . We collect the parameters into a vector A 
with components A> where j = 1..N and N is the total number of parameters. From Bayes' 
theorem 

where p(M\I) is the prior probability density, p(D\M, I) is the likelihood function, and 
p(D\I) = J p(M\D, I) dA is a normalization factor. Our choice of priors is described in the 
next section. 

MCMC is an efficient means of calculating p(M\D, I) whereby one constructs a sequence 
or "chain" of models whose density in parameter space is proportional to the posterior PDF 
provided the chain is long enough to have fully explored all "important" regions of parameter 
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space. Marginalization - that is, integration over a subset of parameters - is trivial; simply 
project the chain onto the appropriate subspace and compute the density of points. Likewise, 
the PDF of any model-dependent quantity is obtained by making a histogram of the quantity 
over the chain of models. 



Our Markov cha in is constructed via the Metrop olis-Hastings algorithm ([Metropolis et al 



19531 ; lHastingsl 1 1970 ) as outlined in [Gregory] (120051 ). The chain of models is described by 



the sequence Aj, i — 0, 1,2... . We begin with a starting point A . A candidate for A 1 is 
chosen according to the jumping rule (also known as the proposal distribution), g(A!|A ). 
Ai is accepted with probability equal to min{l, r} where 

(25) 



p(Ai 


DJ)q(A 


Ax) 


P(A 


D,I)q(A 1 


Ac) 



Otherwise, A 1 is set equal to A . The process is then repeated for A 2 . 

The success of an MCMC analysis rests, by and large, on choosing a suitable jumping 
rule. If the step size from A n to a candidate for A n+ i is too small, the chain will move slowly 
through parameter space. On the other hand, if the step size is too large, most attempts 
to find a new point in parameter space will fail. In either case, exploration of parameter 
space, often referred to as mixing, can require a prohibitively large amount of computing 
resources. Ideally, the jumping rule is shaped like the PDF but scaled by a factor 2A/\/~N 



( IGelman. Carlin. Stern, fc Rubin! Il995l ) which explains why it is often referred to as the 



proposal distribution. 

In this work, we take q to be a multivariate Gaussian so that 

A n+1 = A n + D G (26) 

where G is a vector of length TV whose components are Gaussian random variables with unit 
variance and D is a user-specified N x N matrix. Since neither D nor G depend on the 
model parameters, q (A n+1 |A n ) = q (A n | A n+1 ) and therefore r = p (A n+1 \D, I) jp (A n |Z?, I) 

We begin with a simple ansatz for the proposal distribution in which D is a diagonal 
matrix whose non-zero elements are given by our best guess for the variance of each of the 
model parameters multiplied by 2A/y/N. From a short chain of a few thousand models we 
estimate the covariance matrix 

i3 ~ (A*)(A>) [ ' 

where (...) denotes an average along the chain. Our improved expression for the proposal 
distribution is given by equationl26l with D 2 = (2A 2 /N) C. 
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Each "data point" carries with it a quoted error. Of course, the observer may have 
underestimated the error or there may be features in the data which cannot be explained 
by the model. Both situations can be handled by introducing an u nknown error p arameter 
for each data set which is added in quadrature to the quoted error ([Gregory! 120051 ). For the 
purpose of the MCMC calculation, these error parameters are simply incorporated into an 
expanded definition of A, that is, treated as model parameters. 



5. PRELIMINARIES AND PRIORS 



Simple arguments, based on general features of the Galaxy, provide preliminary esti- 
mates for the model parameters which in turn guide our choices of the prior probabilities 
used in the MCMC analysis. We assume that the priors for each of the model parameters 
are non-zero over a range somewhat larger than the range suggested by these estimates. 
For parameters that correspond to a physical scale (e.g., R e , M d , Vh) we assume a Jeffrey 
prior, essentially, equal probability in logarithmic intervals over the prescribed range. For 
dimensionl ess parameters, such as the halo cusp strength and Sersic index, we assume a flat 
prior. (See iGregoryl (120051 ) for a discussion.) 



The projected velocity dispersion profile to ward the Galactic bu lge, cr p (R), reaches a 
peak value of ~ 130 km s -1 at a radius ~ 200 pc (ITremaine et al.ll2002l ). On the other hand, 
esti mates of the hal f -light or effective radius of the bulge, R e , range from 570 pc to 920 pc 



sec 



Tremaine et al.l (120021 ) and references therein). 



The projected velocity dispersion of the iPrugniel fc Simienl (119971 ) profile exhibits a 
similar structure to that of the Milky Way: <J P (R) is non-zero at R = 0, r ises to a peak value 
of apea k at a radius -R pea k and th en decre a ses w ith ra dius (s e e Fig ure 10 of IPrugniel &: Simien 
(119971 ) as well as earlier work by iBinneyl (Il980l ) and ICiottil (Il99ll )). The dimensionless ratios 
Rpea,k/R e , o'peak/c'b and M bulge /of R e are functions of n as shown in Figured! From the figure 
we deduce that for the Milky Way, n is less than 2, R e is between 0.57 kpc and 0.92 kpc, 
is between 340 km s -1 and 400 km s -1 , and M b is between 1 x 10 10 M - 3.4 x 10 10 M . 



Binney. Gerhard, fc Spergell (119971 ) constructed a model for the luminosity density of 



the Galaxy to fit data from the DIRBE experiment. Their model consisted of a triaxial 
bulge and doubl e expo nential disk with bulge-to-total luminosity ratio of 0.16. Subsequently, 
Malhotra. et al.l (119951 ) derived a total L-band luminosity for the Milky Way of 7.1 x 10 10 L Q 
with 1.1 x 10 10 L Q attributed to the bulge. 

Stellar population s ynthesis models provid e estimates for the stellar mass-to-light ratios 
in different wavebands (IBell fc de Jong! (120011 ) and references therein). The L-band stellar 
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mass -to-light ratio for the disk is expecte d to be in the ra n ge 0.5 — 0.65 in solar units (jde Jong 
20071 ) assuming a scaled Salpeter IMF (jBell fc de Jongi (120011 )) and the Pegase population 
synthesis model. The mass-to-light ratio for the bulge could be somewhat higher. On the 
other hand, since our model does not include a separate gas disk, the effective mass-to-light 
ratio for the disk must account for any gas and should therefore be high er than the value for 
a pure stellar disk. The local stellar-to-gas ratio is ~ 1.6 (see Table 1 of iBinney fc Tremaine 
(119871 )) and therefore the effective (M/L) d might be closer to 1. Together with our estimate 
for the disk luminosity, we conclude that Md is between 3 x 10 10 and 6 x 10 10 M & . 



Reidl (119931 ) reviewed estimates of the distan ce from the Sun to the G alactic center and 
concluded that Rq = 8.0 ± 0.5 kpc. More recently lEisenhauer et al.l (120031 ) observed the star 
S2 in orbit about the Galaxy's massive central black hole using the ESO VLT and found 
R = 7.94 ±0.42 kpc. 



Sackettl (119971 ) reviewed estimates of the radial scale length of the Galactic disk and 
found Rd = 3.0 ± 1 kpc. More recent estimates show a similarly large spread in values. 
Zheng et al.l (120011) fou n d Rd = 2.75 ± 0.3 kpc from HST observations of M dwarfs while 
Lopez- Corredoira et al.l (120021 ) found Rd = 3.3^4 kpc from an analys is of old stellar pop- 



ulations using 2MASS survey data. As emphasized by ISackettl (119971 ). the ratio Ro/Rd is 
observationally more secure than Rd- The estimates collected in her review show Ro/Rd 
between 2.7 and 3.5. 

The disk scale height parameter, hd, is more difficult to constrain since the Galactic disk 
comprises at least three distinct components, the gas disk, the thin disk, and the thick disk 
whereas our model has a single-component disk. Multi-component disks will be incorporated 
into future versions of the model but for the ti me being, hd m ust represent the vertical mass 
distribution of all disk-like components. From ISackettl (119971 ) we surmise that hd is between 
0.2 and 1 kpc. 



The radial velocity dispersion in the solar neighborhood is 36±5.4 km s 1 (IBinney fc Merrifield 



19981 ). Together with estimates of the radial scale length of the disk and with equation [T2| 
this translates into a range of possible values for ctro- 

We allow the halo parameters to vary over a wide range of values. For example, we as- 
sume that the prior probability distribution of 7 is uniform between and 1.5 and non-zero 
otherwise. An alternative approach is to use cosmological models of halo and galaxy for- 
matio n to guide ones choice of the halo parameters (see, for example, IValenzuela fc Klypin 
(120031 )). but given uncertainties in the exact nature of adiabatic compression, variations 
among halo profiles found in different simulations, and possible discrepancies between halo 
profiles as inferred from observations and those found in simulations, we take a more conser- 
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vative approach. Furthermore, since the data do not probe the Galaxy beyond 100 kpc we 
set r h = 100 kpc and 5r h = 5 kpc though no physical meaning should be ascribed th these 
values; can be increased without changing the fit to the data. 

Our choices for the prior probabilities of the model parameters are found in Table 1. 
In addition to the parameters listed in Table 1, we include unknown error parameters (see 
previous section) for the terminal velocity, outer rotation curve, bulge dispersion, and surface 
brightness profile data sets. Thus, our MCMC analysis is run with seventeen parameters. 



6. MCMC RESULTS 

We generate five Markov chains of length 2 x 10 4 . Each chain has roughly 3400 distinct 
members corresponding to an acceptance rate of 17%. We can calculate the PDFs for 
different quantities using the five chains individually and in combination. If mixing has been 
achieved then the results will be the same to within statistical uncertainties. In Figure [5] we 
illustrate that this is indeed the case by plotting the average values with la error-bars for the 
a selection of six model parameters calculated for each of the five chains. For convenience, 
the values are normalized by dividing by the overall average. 



6.1. Selected Models 



In Figure [3] we show the rotation curve and density profile for a model from our MCMC 
series with n ~ 1 and 7 ~ 1. Also shown are the target bulge and halo profiles (equations [7] 
and [9]). We see that the model profiles are very close to the target ones. By comparison the 
halo profiles in the Widrow & Du binski models differ s ignific antly from the target (NFW) 
profile (see Figure 2, top panel, of IWidrow fc Dubinski! (120051 )). 



We display results for the surface brightness profile (Figure 4), the terminal velocity 
(Figure 5), and the bulge line-of-sight velocity dispersion (Figure 6) for three models from 
our MCMC analysis. We choose models with 7 ~ 1 and Sersic parameters n ~ 0.6, 1 and 2. 
We also include a model with n = 4, that is, with essentially an i? 1//4 -law bulge. Since our 
MCMC run does not find any models with a Sersic parameter this large we generate this 
model by fixing n = 4 and allowing the remaining parameters to vary. 

While suitable models are found for n between 0.6 and 2 this is not the case for 

n = 4. This result is in agreement with studies of bulges in late-type spiral galaxies 

(|Andredakic. Peletier. fc Balcellall995l : ICourteau. de Jong, fc Broeilall996l ; iMacArthur. Courteau. fc Holtz 
20031 ) and suggests that the Galaxy has a pseudo-bulge rather than a classical bulge (see 
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Kormendy &: Kennicuttl (120041 ) for a review) 



The models clearly have a difficulty reproducing the shape of the line-of-sight velocity 
dispersion profile. In particular, the local velocity minimum found in the data at R ~ 
3pc is much deeper than is allowed by the models. The discrepancy may indicate that 
the density profile of the Galactic bulge is significantly different from the one proposed by 



Prugniel fc Simienl (119971 ) or that velocity-space anisotropy and deviations from spherical 



symmetry are necessary to model the dispersion profile in the in nermost region of the bulge 
jSpergel. Malhotra. & Blit3ll996l : lFreudenreichlll998h . 



6.2. Statistical Overview 

The maximum a posterior values and boundaries of the 68.3% credibility regions for the 
input parameters and calculated quantities are given in Tables 2 and 3 respectively. Not 
surprisingly, and 7 exhibit the largest fractional uncertainties. Most of the data used in 
this study pertains to the region of the Galaxy near and inside the Sun's orbit about the 
Galactic center; Equation [21] provides a rather weak constraint on the mass at large radii 
(or equivalently, the circular speed at these radii). We note that ah and 7 are correlated in 
the sense that models with larger values of ah have 7 closer to 1 — large const ant- density 
cores (7 ~ and ~ 20 — 30 kpc are disfavored by the data. 

In Figure 7 we show the PDFs for the disk and bulge masses as well as the halo mass 
within 10 kpc and 10 kpc. We also show the pseudo-likelihood function for the twenty-two 



models considered in lDehnen fc Binneyl (119981 ). Specifically, we plot 



C DB = e -(*DB-xLn) (28) 
where Xdb * s from their Table 3 with Xmm se ^ to the value for their b est-fit model. Our 



result s for the disk, bulge, and halo masses are consistent with those of iDehnen &: Binney 



( 119981 ). Figure [7] illustrates the advantages of the MCMC method. The Dehnen & Binney 
analysis involves twenty-two separate maximum-likelihood calculations characterized by the 
authors' ad hoc choices of fixed and free parameters. MCMC, on the other hand, yields 
the full multi-dimensional posterior probability function from which PDFs for one or more 
parameters are easily obtained. 

Figure M provides a contour plot of the PDF in the R — Rd plane. Our results are 
consistent with previous estimates of these two quantities. The plot also shows the general 
trend that models with higher values of Rd tend to favor higher values of R$. 



In Figure 9 we show the PDFs for n and 7. As noted above, our analysis clearly favours 
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bulges with surface brightness profiles closer to an exponential than to de Vaucouleurs i? 1 / 4 - 
law. Our results allow for a dark halo with a wide range of inner logarithmic density slopes 
which includes the NFW value as well as steeper and shallower indices. 



6.3. Comparison with Published Milky Way Models 

In Figure [101 we compare values for the disk and bulge masses from our Markov chain 
analysis with those from a num ber of popular Milky Way models . One of the earliest 
mass models was constructed by iBahcall. Schmidt, fc Soneiral (119831). While they focus o n 
constraining the Galactic spheroid through star counts (see also IBahcall fc Soneiral (119801 )) 
they also fit the Galactic rotation curve by modeling the mass distribution of the disk, 
bulge, and dark halo. Their choice for the disk and bulge masses (Md = 5.6 x 10 10 M & and 
Mb = 1.1 x 10 10 M Q ) is represented in Figure [TUl by the solid triangle. 



From Shuttle IRT observations, iKent. Dame, fc Faziol (119911 ) derived total K-band lu- 
minosities for the disk and bulge of Ld = 4.9 x 10 10 L Q and Lb = 1.1 x 10 10 L , respectively 
(Note that these values differ slightly from the values qu oted i n the o riginal paper because 
a different value for the K magnitude of the Sun is used. iKentl (119921 ) went on to construct 
disk-bulge-halo mass models based on these results and was able to fit the Galactic rotation 
curve for three choices of the disk mass-to-light ratio: a maximal disk model ((M/L) d = 1.3), 
a high disk model ((M/L) d = 1.0), and a low disk model ((M/L) d = 0.68). In each case, the 
bulge mass-to-light ratio was (M/L) b = 1.0. Kent's models are shown in Figure [TOl as blue, 
red, and green stars. 



Klypin. Zhao, fc Somervilld (120021 ) constructed models for the Milky Way and An- 
dromeda galaxies based on cosmologically motivated theories of disk formation. Their fa- 
vored model for the Milky Way has M d = 4.0 x 10 10 M and M b = 0.8 x 1O 1O M and is 
represented in Figure [10] by the solid square. 



Since its discovery (llbata. Gilmore. &: Irwinlll994l . Il995l ). the Sagittarius dwarf galaxy 



has held the promise of providing useful constraints on the size and shape of the Milky Way's 
dark halo. The usual approach is to simulate the tidal disruption of Sagittarius as it passes 



and kinematic observations of the observed tidal streams 



(see 


Law, Johnston, & Maiewski 


Johnston et al. 


(1999 


) introduced 



a m odel for the Gal a ctic p otential which has now become a standard for work in this area. 
The Ijohnston et al.l Jl999h values, M d = 1.0 x 10 11 M Q and M b = 3.4 x 10 10 M , are shown 



in Figure [TD] as an open square. 
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Our models occupy a smaller region of Md — Mb parameter space than is spanned by 
popular models from the literature. The iJohnston et all (119991 ) choices for Md and Mj, 
are inconsistent with our results by factors of 2.5 and 4, respectively. More to the point, 
their choices fa ll well outside the region of accep table mode ls. The choices for disk and 



bulge masses in iBahcall. Schmidt, fc Soneiral (119831) andlKent Jl992h are consistent with our 



results as is the preferred model from Klypin. Zhao, fc Somervillel ( 2002 ). 



6.4. Implications for Dark Matter Detection Experiments 



Models of the Milky Way are an essential ingredient in the analysis of dark matter 
detection experiments. For example, microlensing experiments, which attempt to measure 
the distribution of massive compact halo objects (MACHOs) along various lines-of-sight 
through the halo, require a model for the MACHO component of the halo as well as the 
distribution of stars in the disk. 

In analyzing their 5.7 year data set, the MACHO experiment considered a wide range of 
Galactic models (lAlcock et al.ll2000l ). Here we focus on their "standard model". This model 
includes both a thin and thick disk, each with radial scale length Rd = 4kpc and total disk 
mass is Md = 4.5 x 10 10 M Q . The halo is modeled as a cored isothermal sphere with a density 
profile given by 



P{r) = Po 



a' 



+ RI 



a 2 + r 2 



(29) 



where, for their standard model, lAlcock et all (120001 ) assumed Rq = 8.5 kpc, a = 5kpc 
and po = 0.0079 M Q pc -3 . Though they did not model the bulge explicitly, we can infer its 
mass through the requirement that the circular rotation speed of the Galaxy at r = R Q is 
~ 220kms~ . Doing so yields Mb ~ 2.6 x 1O 1O M for their standard model (open triangle 
in Figure [TO]) . Evidently, the bulge mass is inconsistent with our results by a factor of 3. 

Terrestrial dark matter detection experiments aim to observe elementary particle dark 
matter candidates (e.g., WIMPs or axions) as they interact with a detector on Earth. These 
experiments are therefore sensitive to the local density and velocity distribution of dark mat- 
ter particles. Estimates fo r the local dark matter density ran g e from 0.005 — 0.02 Mnpc "' 3 
0.2- 0.8Ge Vcm- 3N ) (See bahcaU. Schmidt, fc Soneiral Jl983h : ICaldwell fc Ostrikel 



( 



1981) 



?urnerl ( 119861 ) ; IBergstrom. Ullio. Sz Buckley! (119981 ) ) . The range quoted above is from lBergstrom. Ullio. fc B 
19981 ) where a variety of halo profiles were considered. 



We find (Table 3) pi oca i = 0.0080 ± 0.0014 M pc" 3 and ai oca i = 240 ± 23kms" 1 . The 
PDFs for these two quantities are plotted in Figure [H Our values for the mean and lower 
bound for pi oca i are consistent with the values found in IBergstrom. Ullio. fc Buckley! (119981 ) 
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though our analysis suggests that their upper bound is too high by a factor of two. Our mean 
value for cr local is lower than the standard value by about 30 km s -1 though the standard value 
still falls within the range of acceptable models. Our models tend to favour values for M wo 
at the low end of the range found in equation [5TJ Inspection of a scatter plot of models in 
the Mioo — cxiocai plane reveals a clear correlation between the two quantities — models with 
values of Mi o closer to the central in equation [21] have values of cri OC ai closer to the standard 
value. 



6.5. Connection with Cosmology 



Klypin. Zhao, fc Somervilld (120021 ) construct models for the Milky Way based on stan- 
dard galaxy formation theory. In particular, they use model halos based on predictions 
from cosmological simulations with the further assumption that the halos undergo adiabatic 
contraction in response to the formation of the disk and bulge. 

Though our MCMC analysis does not include cosmological constraints we can test 
whether our models are consistent with the standard cosmological paradigm a posteriori. To 
this end, we calculate the virial radius, -R v ir, virial mass, M v i r , and concentration, c v ; r , for all 
of the models in our MCMC series. By definition, the mean density inside i? v ; r is a factor 
A, 



greater than the background density, p m . That is, 

M vir = M (-Rvir) = -^XwP m Rt 



4vr 
~3~ 



(30) 



where M vir is the mass interior to i? vir . A vir depends on the cosmological model; In what 



Pm/P, 



follows, we assume JX^ 
H = 70kms _1 Mpc (Tegmark et al 



1998h (see also iBullock et all (120011 ) and 



0.3 where p cr i t = 3H 2 /8n is the cr itical density and 



2004 



i. With these values, A v j r ~ 340 (iBryan &: Norman 
Wechsler et al.l (120011 )). The concentration, c vir , is 



defined as the ratio of R vir to the halo scale length R s , the latter identified as the radius at 
which the logarithmic slope of the halo density profile equals -2. For a halo profile given by 
equation @ ]? s = (2 - 7) a^. 

In Figure 12, we show the PDF for our MCMC series projected onto the M vir — c v - ir and 
-Rvir — Cvir pl anes. Also shown are the low-con centration, high-concentration, and favored 
models from iKlypin. Zhao, fc Somervilld (120021 ) . In the upper panels, i? v i r , M vir , and Cvir 
are calculated for the actual halos used in our models. The i mplication would seem to be 
that o ur models are more concentrated than the ones assumed in lKlypin. Zhao, fc Somerville 
J2OO2I). 



Recall, however, that IKlypin. Zhao. &: Somervilld (120021 ) incorporate adiabatic contrac- 
tion into their models. In order to compare our halos with theirs, we should assume 
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they too have undergone adiabatic contraction. We should therefore adiabatically "de- 
contract" our halos and then calculate Cy ir . We have done this using the usual assumptions 
( IBlumenthal et al.lll986l ; iFlores et al.lll993l ); The system is treated as being spherically sym- 
metric. Initially, the baryons and dark matter are well-mixed. The baryons cool and form a 
disk and bulge while the halo particles respond adiabatically to the changing gravitational 
potential. Moreover, halo particles are assumed to follow circular orbits which do not cross 
as their orbits shrink. Under these assumptions, the quantity rM(r) remains constant and 
one can calculate the initial radii of the dark matter particles given the final structure of the 
disk, bulge, and halo. 

The results of this analysis are shown in th e lower panels of Figure 12. Our m odels are 
now in excellent agreement with those found in iKlypin. Zhao, fc Somervilld (120021 ) . 



Alam. Bullock, fc Weinberg! (120021 ) proposed Ay/ 2 , the mean density within the ra- 
dius Rv/2 in units of p cr j t , as a measure of the halo core densities. Rv/2 is defined as 
the radius at which the rotation curve reaches half its maximum value V max . The quantity 
(Ay/ 2 /87r 2 ) is equal to the numb er of rotation periods per Hubble time at the radius R 



V/2- 

Alam. Bullock, fc Weinberg! (120021 ) found values of Ay/ 2 m the range 5 x 10 4 — 5 x 10 6 for 
dark matter-dominated galaxies, considerable scatter and no correlation between Ay/2 and 

^max- 

The PDF for our MCMC series in the Ay/ 2 — V mSuX plane is s hown in Figure 13. The range 
of valu es for Ay/ 2 is certainly consistent with those found in the lAlam. Bullock, fc Weinberg 
(120021 ) survey though some models have higher values — possibly reflecting that influence 



of baryons on the Milky Way's dark halo. The considerable range indicates that the halo 
concentration is poorly constrained by the data used in our analysis. 



7. BAR FORMATION 



Near-IR photometry, gas and stellar kinematic measurements, and observations of grav- 
itational microle nsing ev e nts a ll stro ngly suggest th at the Milky Way is a barred galaxy 
(see reviews by iKuijkenl ( 119961 ) and iGerhardl (119961 )). A bar represents a strong depar- 
ture from axisymmetry and adds considerably to the challenge of modeling the Galaxy. A 
promising avenue is to use N-body simulations to follow the development of bar s and spiral 
structure in an initially axisymmetric, equilibrium model. lOstriker fc Peebles! (119731 ) and 
Sellwoodl (119851 ) provide early examples of this approach. They were interested in stabilizing 
their Galactic models to avoid bar formation, the former proposin g an unseen dark halo 
and the latter stressing the importance of the bulge. More recently, iFuxl (119971 ) attempted 
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to construct self-consistent models for the Milky Way's bar by evolving unstable axisym- 
metric models and comparing the results with observations from the DIRBE experiment 
jDwek. et alJll995h . 



Bars, at least in ideali zed, initially-axisymmetric disk galaxies, form through swing 
amplification (jToomrelll98ll ). (Whether bars in real galaxies form through this process or 
through some more complicated process during the formation of the galaxy itself is another 
matter.) In swing amplification, leading waves propagate outward and are amplified into 
trailing waves (and contained by the outer Lindblad resonance). Trailing waves then wind 
up. Within linear theory, if there is an inner Lindblad resonance (ILR), it absorbs the tightly 
wound trailing waves and thereby acts as a barrier preventing further growth. In the absence 
of an ILR, the trailing waves propagate through the center of the system and transform into 
leading waves. This feedback loop can lead to growth of a bar-like perturbation. 



The Toomre Q ( )Toomrelll964l ) and Goldreich-Tremaine X ( iGoldreich fc Tremaindll978 



19791 ) parameters are the two most widely used diagnostic quantities for studying a galaxy's 
resistance to the bar instability. Q measures the kinetic "temperature" of the disk; stellar 
disks with Q < 1 are unstable to local gravitational instabilities. X indicates a disk's 
susceptibility to global instabilities. (X depends on the azimuthal mode number, m, of 
the perturbation. Here, we take m = 2 since we are interested in bars.) For disks with 
X < 3, the gain of the swing amplifier is large and bars are more likely to form. In general, 
the greater the contribution of the disk to the gravitational force felt by disk particles, the 
smaller the value of X. 

Higher values of Q and X make a galaxy more resistant to the bar instability. The 
existence of an ILR barri er from a dense bulge or cuspy halo can also prevent the instability 
( jSellwood fc Evans! 120011 ). However, a bar can still form even if the galaxy initially has an 
ILR barrier. For example, interactions between halo substructure and the bulge or cusp may 
temporarily lower the ILR barrier an d trigger bar formation in an otherwise stable galaxy 
( IGauthier. Dubinski. fc Widrowl 120061 ) . Non-axisymmetric structure in the disk might also 
jostle the cusp and remove the ILR barrier, if only temporarily. In short, the notion of 
an ILR barrier assumes linear perturbation theory; non-linear disturbances may be able to 
overcome or disrupt the barrier and initiate the bar instability. 

Our models provide a natural starting point for investigations of bar formation. The 
models generated by our MCMC run all yield acceptable fits to the observational data 
yet have very different stability properties. Figure 14 shows a contour plot of the model 
distribution in the Q — X plane. (Of course, both Q and X are functions of radius. Here we 
use their minimum values.) Also shown is the distribution of models in the X' — X plane 
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where 

xrl _ W total("^) (n-i\ 
^disk W R=2.2R H 



id 



is a measure of the disk's contribution to the gravitational force necessary to keep a particle in 
a circular orbit at a given radius. Here v to tai is the circul ar rotation speed at cylindrica l radiu s 



R and t>di S k is the contribution to t> to tai from the disk (IDebattista fc Sellwoodlll998l . 120001 ) . 
The radius R = 2.2R& is where Udisk reaches a maximum, assuming an exponential disk. The 
tight correlation between X and X' indicates these two quantities are interchangeable for 
most purposes. 

We consider twenty-five models from our MCMC run which cover most of the area within 
the 95% likelihood contour in the Q — X plane. The models in this study are depicted as 
dots in Figure [TH while their circular-speed curves are shown in Figure [15j In all cases, the 
circular speed reaches a peak value of approximately 220 km s -1 at a radius between 5 — 8 kpc 
and then declines slowly. 

Figure [T6l shows the behavior of the functions Q and Q ± k/2 where Q is the angular 
velocity and k is the epicyclic radial frequency. These functions should be compared with 
the pattern speed, Q p , of an emerging bar or spiral density wave. In particular, Q = Q p 
indicates the position of co- rotation while Q p = Q — k/2 indicates the position of the ILR. 
As we will see shortly, all of the bars that form in our simulations have initial pattern speeds 
in the range of ~ 40 — 50 kms" 1 kpc" 1 . Thus, about two thirds of the models in this study 
have initial ILRs. Note that in these cases, the rise in the Q — ft/2-curve is very rapid and 
the radius of the ILR is only a few hundred parsecs. 

Each of the twe nty-five models depicted in Figure [TH are evolved for 5 Gyr using a 



parallelized treecode (IDubinskil Il996l ). Simulations have 800K disk, 200K bulge, and 1M 
halo particles. The particle Plummer softening length is e = 25 pc and the simulations are 
run for 10 4 equal timesteps of length At = 0.5 Myr. We generate surface density maps for 
the face on view of every tenth timestep and determine the strength and pattern speed of the 
bars that form. An animation depicting the evolution of the grid of simulations is available 
at |http: / /www.cita.utoronto.ca/ ~dubinski/DynamicalBlueprints/| 

In Figure [T7| we plot the amplitude of the bar-strength parameter 

1 



A 2 = 
2 N 



E 



e 2i<t>j 



(32) 



as a function of time. Here <pj is the usual azimuthal angle of the j th disk particle and iV is 
the total number of disk particles. All of the models form bars, with the possible exception 
of the Q = 2, X = 4.5 model (upper right-hand corner of the figures). Bars form almost 
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immediately in models with low values of Q and X. Indeed, the models in the lower left- 
hand region of the figures are violently unstable and radically transformed by both local and 
global instabilities. On the other hand, bars can take several billion years to form in models 
with higher values of Q and especially higher values of X. 

Interestingly enough, the growth-rate of bars does not vary smoothly across the grid of 
models. In particular, models in the central column of the grid (Q ~ 1.5) with higher values 
of X quickly develop bars. In general, these models have less massive and less dense bulges 
and no ILRs (see Figures 15 and 16). Hence, they are immediately vulnerable to Toomre's 
swing amplification mechanism. 

Virtually all of the models with ILRs eventually form bars though the onset of the 
instability is typically delayed. Two effects, one numerical and one physical, can lead to bar 
formation in the presense of an initial ILR. First, two-body relaxation may be important at 
radii near the ILR even with 10 6 halo particles. (Recall that in some models, the radius of 
the ILR can be as small as 100 pc.) As the central halo relaxes, the peak of the Q — ft/2-curve 
can drop below the initial pattern speed. If this occurs, the ILR barrier disappears and the 
bar forms. 

A second, more physical, explanation is that swing amplified spiral waves disturb the 
central region of the galaxy causing the ILR barrier to disappear, at least temporarily. Other 
non-linear disturbances, such as interactions between the disk and halo substructure, may 
also cause the ILR barrier to disappear. Simply put, the theory of the ILR barrier may not 
apply to situations where non-linear perturbations a re present. We n ote that the importance 



of nonlinear effects for bar formation is discussed in ISellwoodl (119891 ). 



The evolution of the bar length, R^, and pattern speed, £l p , are shown in Figures 
and dHJ To estimate the bar length, we first construct the axis ratio profile of the isodensity 
contours for the disk. The axis ratio profile is defined to be the axis ratio, b/a, as a function 
of the semi-majar axis, a. In general, b/a goes through a minimum with a typical value of 
{b/a) mm — 0.4 before rising abruptly in the transition from the end of the bar to the outer 
disk. As a heuristic measure of Rf, we use the major axis length of the isodensity contour 
beyond the 6/a-minimum and where b/a ~ 0.6. The pattern speed is given by the time- 
derivative of the phase angle 0o = tan _1 [Im(y4 2 )/i?e(A2)]/2. We insure that 0o is sampled 
with enough time-resolution to avoid aliasing and to obtain a reasonable estimate of the 
pattern speed from A<p/At. In Figure fl~9| we show the pattern speed once the bar is easily 
detected above the noise. 

Though the bars in our survey appear at different times, their initial pattern speeds are 
always in the range of 40 — 50kms _1 kpc -1 . Angular momentum transfer to the halo causes 
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fl p to decay to about 20 km s _1 kpc~ x over a few Gyr with the decay rate being slightly larger 
in models with higher mass halos. As expected, the bar lengths are generally less than the 
co-rotation radius though transient long bars with R^ « 10 kpc form in models with small Q 
and X before collapsing on themselves. For most models, bars grow monotonically in length 
as the pattern speed declines. 

Bar formation alters the structure of the model and it is therefore natural to ask whether 
the evolved models still satisfy the original observational constraints. Though we will leave 
the details of such an analysis for future work, we include, in Figure [201 the evolution of 
the surface density profiles for our twenty-five simulated systems. We see that the sur- 
face density profiles of models with low values of Q and X are dramatically deformed; the 
redistribution of mass is so violent that the models almost certainly do not satisfy our ob- 
servational constraints. On the other hand, models where the bar forms relatively late in 
the simulation show little evolution of the surface density profile. Their structure, at least 
in an azimuthally-averaged sense, remains largely unchanged. 

The Galaxy is known to h ave a bar and estimates of it s lengt h and pattern speed can 



be compared with our results. iBinney. Gerhard, fc Spergell (119971 ) find Rb = 3.5 k p c and 



Qp ~ 60 — 70kms 1 kpc 1 in their analysis of DIRBE photometry while iDehnenl ( 19991 ) 



find a similar bar length but lower pattern spe ed (Q p = 53 ± 3kms *kp c x ) using the 



velocity distribution of solar neighborhood stars. I Weiner &: Sellwoodl (119991 ) model the gas 
kinematics of the inner galaxy and find Q p ~ 42 km s _1 kpc -1 . The fact that nearly all of the 
models in our study are bar unstable and have initial pattern speeds near the range of the 
inferred values is promising. We note that certain models can be excluded such as those with 
very small values of Q and X. If we assume that the Galaxy's bar has formed very recently 
than many of the models have the correct combination of bar length and pattern speed to 
match the observations. Furthermore, for larger Q and X, the change in the disk's radial 
profile in response to the bar is small. These models may well provide a good barred model 
of the Galaxy. It should be noted that in all of these models, the pattern speed declines to 
20 — 30kms _1 kpc -1 within a few Gyr after the bar forms. If we take these models seriously 
as reasonable facsimiles of the Galaxy, then we must conclude that Galactic bar formed within 
the last 1-2 Gyr. 



8. CONCLUSIONS 



We have introduced a dynamical model for late-type galaxies that incorporates our 
current understanding of disk-bulge-halo systems. In particular, the bulge has a Sersic 
surface density profile and the halo has a central density cusp. 
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We have carried out an MCMC analysis of dynamical models for the Milky Way using 
a variety of kinematic and photometric constraints. The results are presented in the form 
of PDFs for both input parameters and derived quantities. The MCMC analysis provides 
a picture of the distribution of models in parameter space that is more complete than can 
be obtained by other approaches. Avoided is the awkward procedure of fixing a subset of 
parameters while allowing the remaining parameters to vary in some minimization scheme. 
Instead, a sequence of models is generated which contains all of the desired information. 
Marginalization over a subset of parameters is accomplished by simply projecting the model 
distribution onto the appropriate parameter subspace. 

Our analysis suggests that the Milky Way has a pseudo-bulge with a Sersic index of 
1.3 ± 0.3. Our results for the masses of the disk, bulge, and halo are consistent with those of 



Dehnen fc Binneyl (119981 ) but call into question choices for these qua ntities in some popular 



models from the literature. For example, the disk and bulge masses in I Johnston et al.l (119991 ) 



are entirely inconsistent with our res ults. The inferred b ulge mass for the standard model 



used by the MACHO collaboration (lAlcock et al.l 120001 ) is inconsistent with our findings 



by a factor of 2.5. On the other hand, the standard values used by terrestrial dark matter 
detection experiments for the local dark matter density and velocity dispersion are consistent 
with our results. 

A weak point of our analysis is the inability to tightly constrain the halo mass at large 
radii. Planetary nebulae, globular clusters, and satellite galaxies may be used as tracers of 
the Galactic potential. Dynamical models for the tracer populations are required to properly 
model kinematic data. In principle, it is straightforward to construct such models but there 
are subtle issues. Previous studies showed that velocity anisotropy in a tracer population can 
affect interpretation of kinematic data. Velocity anisotropy requires a DF that depends on at 
least one integral of motion in addition to the energy. The standard practice is to use the total 
angular momentum, L, but since our models include a disk, L is not conserved. (Previous 
analyses side-stepped this issue by using a spherically symmetric Galactic potential.) 

Our MCMC analysis provides an ideal starting point for studies of disk stability and 
bar formation in that we have some 10 5 models, each of which can serve as initial conditions 
for a numerical experiment. We have performed a suite of simulations which focuses on the 
susceptibility of the disk to bar formation as a function of the stability parameters Q and 
X. Many of the models provide a good match to the inferred properties of the Galactic bar 
with the proviso that the bar has formed recently. 
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Fig. 1. — Cha r acteri stics of the line-of-sight velocity dispersion profile for the 



Prugniel fc Simienl (119971 ) density profile as a function of Sersic index n. The curve a p (R) 



rises to a peak value of f pea k at a projected radius -R pea k- Shown are the dimensionless ratios 
Rpeak/Re (upper panel), fpeak/^b (middle pa nel), and GMhu^ p. /o'hR e (lower panel) as func- 
tions of n. Data for the Galaxy compiled by iTremaine et al.l (120021 ) suggest that R pea k/Re 
lies in the range 0.22 — 0.36 (horizontal lines). Vertical lines then delineate range of values 
for n. 
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Fig. 2. — Average values and la error bars for a selection of six model parameters as 
calculated from each of the five separate Markov chains. Values of the parameters are 
normalized by dividing by the overall average. 
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Fig. 3. — Rotation curve (top panel) and density profile (bottom) as a function of r for a 
model with n ~ 1 and 7 ~ 1. Red (dashed) line is for the bulge; blue (dotted) line is for 
the disk; green (long-dashed) line is for the halo. In the top panel, the solid (black) line 
shows the total rotation curve. In the lower panel, we plot r 2 p for the bulge and halo and 
h^rp for the disk, quantities proportional to the mass in radial bins. Also shown in the lower 
panel (thin lines) are the "target" density profiles for the bulge and halo (equations [7] and H] 
respectively). 
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Fig. 4. — Surface brightne s s pro file as a function of I for four models and for data from 



Spergel. Malhotra. &: Blitzl (119961 ). Three models are chosen from our main MCMC run: 
(n, 7) ~ (0.6, 1) — red curve; (n, 7) ~ (1, 1) — magenta curve; (n, 7) ~ (2, 1) - - blue 
curve. Also shown (green curve) is a model chosen from a run where n is fixed to the "de 
Vaucouleurs" value, 4 and 7 = 0.86. In the top panel, the thin curves show the separate 
contributions of the disk and bulge. Lower panel shows the residuals between the models 
and the data. 
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Fig. 5 . — Terminal velocity as a function of sin I for four models and for data from iMalhotr, 
(119951 ) Line types and colours are the same as in Figure 4. 
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Fig. 6. — Line-of-sight velocity dispersion toward the bu lge as a function of pro jected radius 
from the Galactic center for models and for data from iTremaine et al.l (120021 ). Line types 
and colours are the same as in Figure 4. 
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Fig. 7. — PDFs for the disk mass, bulge mass and halo mass within 10 and 100 kpc. Solid 
lines show results of the MCMC analysi s . Dot s show likelihood functions for the twenty-two 
models presented by iDehnen &: Binneyl (119981 ) . 
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Fig. 8. — Contour plot of probability distribution function of models in the Rq — Rd plane. 
Solid contours enclose 68% and 95% of the models. Dashed contours enclose 38% and 87% of 
the models. The solid straight (black) line corresponds to Ro/Rd — 3. The dashed straight 
lines correspond to Ro/Rd = 2.7 and 3.3. 
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Fig. 9. — PDFs for the Sersic index, n, and cusp-strength parameter, 7. 
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Fig. 10. — Comparison of our results for the disk and bulge mass with those from pub- 
lished Milky Way models. Cloud of points in the lower right corner are from the MCMC 
analysis. Stars from left to ri ght (b l ue, re d, green) represent, respectively, the low, 
high land maximal models from iKentl (119921 ). The filled triangle represents the model 
from Bahcall. Schmidt, fc Soneiraj ( 1983 ) while the filled square represents the model from 
Klypin. Zhao, fc Somerville! hoO$ \. The open square represents the model advocated by 



Johnston et al. 



(119991 ) in their studies of the tidal disruption of Sagittarius. The open trian- 



gle is the standard model adopted by the MACHO collaboration. 
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Fig. 11. — PDFs for the local dark matter velocity dispersion, o"i OC ai and dark matter density, 
Aocai- Vertical lines indicate the standard values assumed by most terrestrial dark matter 
detection experiments (see text). 
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Fig. 12. — Contour plots of the probability distribution function of models in the M V - 1V — c v ; r 
and R vir — c vir planes. Solid contours enclose 68% and 95% of the models. Stars indicate 
the favored, low-concentrat i on few = 5) and high-concentration (c vir = 17) models from 



Klvpin. Zhao, fc Somervilld (120021 ) . 
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Fig. 13. — Contour plot of the probability distribution function of Ay/2 — V v 
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Fig. 14. — Contour plots of probability distribution function of models in the Q — X and 
X' — X planes. Solid contours enclose 68% and 95% of the models. Dots correspond to 
models used in bar formation study in Section 7. 
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Fig. 15. — Rotation curves for the twenty-five models used in our bar formation simulations. 
Models are arranged so that Q increases to the right and X increases from bottom to top. 
Shown are the total rotation curve (long-dashed line) and contributions to the rotation curve 
from the disk (solid curve), bulge (dotted curve), and halo (dashed curve). 
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Fig. 16. — £1 — nn/m as a function of radius for n = (solid curve), n — 1, m — ±2 (dotted 
curves) . 
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Fig. 17. — Growth of the bar strength parameter, A 2 , as a function of time. 
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Fig. 18. — Bar length, R b as a function of time. 
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Fig. 19. — The bar pattern speed fl p as a function of time. Bars are born with pattern speeds 
fl p ~ 50kms _1 kpc -1 which immediately begin to decay as they transfer angular momentum 
to halos. 
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Fig. 20. — Surface density profile evolution. Three times are shown: the initial time t = 
(black), t = 2.5 Gyr (red), and t = 5.0 Gyr (blue). 
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Table 1. model parameters and prior probabilities 



L_/CArX CAji.J--l.V-' Li\_^X 


prior 


lower hound 


unnor bound 




Jeffreys 


2kms" 1 


6kms _1 


ah 


Jeffreys 


2kpc 


35 kpc 


1 


uniform 





1.5 


M d 


Jeffreys 


2 x 10 10 M Q 


7 x 1O 1O M 


R d 


Jeffreys 


2 kpc 


3.5 kpc 


hd 


Jeffreys 


0.2 kpc 


1 kpc 


°\ro 


Jeffreys 


80 knis" 1 


300 kins" 1 


n 


uniform 


0.6 


2.0 


Vb 


Jeffreys 


150 kins" 1 


400 kms" 1 


Re 


Jeffreys 


0.4 kpc 


1 kpc 


(M/L) d 


uniform 


0.6 


1.2 


(M/L) b 


uniform 


0.6 


1.2 


Ro 


Jeffreys 


7 kpc 


9 kpc 



Table 2. results for input parameters 



parameter average 



ah 

7 

M d 
R d 

h d 
n 

Re 

(M/L) d 
(MjL) h 
Ro 



Note. - Units are: 
kms -1 for velocities, 
10 10 Mq for masses, and 
kpc for distances. 



330±f 5 

0.81 +0 39 39 
a -1+0.53 
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Table 3. results for calculated quantities 



parameter average 
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Note. - Units for p are 
M Q pc~ 3 . Other quantities 
use the same units as in Ta- 
ble 2. 



